%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 遗传算法
% 简介：遗传算法寻找最优路径
% 作者：Zhaojiang
% 日期：2023/10/11
% 企鹅：277746470
% 最短长度：177.2600
% 最短路径：1->2->27->22->25->10->19->16->11->34->13->20->17->24->5->3->29->
%          28->30->7->9->6->4->33->32->23->14->8->21->26->31->18->12->15->1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
return;
%% 数据加载预览
tsp=TSP(importdata('./city_data.mat'));
figure('Name','遗传算法');
tsp.draw_path;
disp 当前路径
tsp.print_path;
pause(2);
close all;
figure('Name','遗传算法','Position',[200 200 1280 480]);
%% 设定算法参数
rng(1234);
elitism = true;% 精英选择
population_size = 200;% 种群大小
generation_size = 5000;% 迭代次数
cross_rate = 0.1;% 交叉概率
mutate_rate = 0.01;% 变异概率
chromosome_size = tsp.city_num;% 计算染色体长度（实数编码）
%% 算法数据初始化
fitness_value = zeros(1,population_size);% 适应度值
fitness_avg = zeros(1,generation_size);% 平均适应度
fitness_sum = zeros(1,population_size);% 总适应度
fitness_best = 1/tsp.path_length;% 当前最佳适应度
best_fitness = zeros(1,generation_size);% 记录最佳适应度
best_generation = 0;% 最佳繁殖代数
population = zeros(population_size,tsp.city_num);
init_path = importdata('./Best_Path_GA.mat');
for i = 1:population_size
    population(i,:) =[1 randperm(tsp.city_num-1)+1];% 随机初始化种群
%     population(i,:) =init_path;% 最优初始化种群
end
population(1:ceil(population_size/100),:) = repmat(init_path,...
                [ceil(population_size/100),1]);
%% 开始迭代求解
for generation=1:generation_size
    % 计算适应度
    % 变量归零
    for i = 1:population_size
        % 适应度计算（根据实际情况调整)
        tsp.update_path(population(i,:));
        fitness_value(i) = 1/tsp.path_length;
        % 保存最佳适应度
        if fitness_value(i) > fitness_best
            fitness_best = fitness_value(i);
            best_generation = generation;
            best_chromosome = population(i,:);
        end
    end
    % 种群排序
    % 冒泡排序总适应度
    for i = 1:population_size
        min_index = i;
        for j = i+1:population_size
            if fitness_value(j) < fitness_value(min_index)
                min_index = j;
            end
        end
        % 如果最小值不是当前值，则交换
        if min_index ~= i  
            temp = fitness_value(i);
            fitness_value(i) =  fitness_value(min_index);
            fitness_value(min_index) = temp;
            
            temp_chromosome = population(i,:);
            population(i,:) =  population(min_index,:);
            population(min_index,:) = temp_chromosome;
        end
    end
    % 计算总适应度
    for i = 1:population_size
        if i==1
            fitness_sum(i) = fitness_value(i);
        else
            fitness_sum(i) = fitness_sum(i-1) + fitness_value(i);
        end
    end
    % 计算平均适应度
    fitness_avg(generation) = fitness_sum(end)/population_size;
    % 记录最佳适应度
    best_fitness(generation) = fitness_best;
    % 逐个选择种群(轮盘赌)
    new_population = zeros(size(population));
    for i=1:population_size
        % 随机选择种群适应度
        select_fitness_sum = rand * fitness_sum(end);
        % 择中法选取最接近随机选择的适应度
        first = 1;
        last = population_size;
        mid = round((first+last)/ 2);
        select_population_id = -1;
        while(first <= last) && (select_population_id == -1)
            if select_fitness_sum > fitness_sum(mid)
                first = mid;
            elseif select_fitness_sum < fitness_sum(mid)
                last = mid;
            else
                select_population_id = mid;
                break;
            end
            mid = round((first+last) / 2);
            if(last-first) == 1
                select_population_id =last;
                break;
            end
        end
        new_population(i,:) = population(select_population_id,:);
    end
    % 更新种群信息
    elitism_population = population(end,:);% 精英种群
    population = new_population;% 选择后的种群
    % 是否保留精英种群
    if elitism
        population(end,:)=elitism_population;
    end
    % 父母交叉染色体(基于位置顺序)
    for i = 1:2:population_size
        fater1_chromosome = population(i,2:end);
        fater2_chromosome = population(i+1,2:end);
        fater2_chromosome_copy = population(i+1,2:end);
        son1_chromosome = fater1_chromosome;
        son2_chromosome = fater2_chromosome;
        cross_position = [];% 子代1交叉位置
        uncross_position = [];% 子代1未交叉位置
        for j = 2:chromosome_size
            if rand < cross_rate
                son1_chromosome(j) = fater1_chromosome(j-1);
                fater2_chromosome(fater2_chromosome==son1_chromosome(j) ) = [];
                cross_position = [cross_position;j];
            else
                son2_chromosome(j) = fater1_chromosome(j-1);
                fater2_chromosome_copy(fater2_chromosome_copy==son2_chromosome(j) ) = [];
                uncross_position = [uncross_position;j];
            end
        end
        assert(length(uncross_position)==length(fater2_chromosome))
        for k = 1:length(uncross_position)
            son1_chromosome(uncross_position(k)) =  fater2_chromosome(k);
        end
        assert(length(cross_position)==length(fater2_chromosome_copy))
        for k = 1:length(cross_position)
            son2_chromosome(cross_position(k)) =  fater2_chromosome_copy(k);
        end
    end
    % 染色体突变（截取片段乱序）
    for i = 1:population_size
        if rand < mutate_rate
            start_position = round(rand*chromosome_size)+2;
            fragment_length = round(rand*(chromosome_size-start_position)) - start_position;
            if fragment_length > 0
                fragment_chromosome = population(i,start_position+1:fragment_length+start_position);
                fragment_chromosome = fragment_chromosome(randperm(fragment_length));
                population(i,start_position+1:fragment_length+start_position) = fragment_chromosome;
            end
        end
    end
    if mod(generation,generation_size/50)==0 || generation ==1
        subplot(1,2,1);hold off;
        tsp.update_path(best_chromosome);
        tsp.draw_path;
        disp(['迭代次数：' num2str(generation) '/' num2str(generation_size)])
        disp(['最短长度:' num2str(1/fitness_best)])
        subplot(1,2,2);hold off;
        plot(1./best_fitness);
        drawnow;
%         pause(0.1);
%         break;
    end
end
%% 处理结果
clf;
subplot(1,2,1);hold on;
tsp.update_path(best_chromosome);
tsp.draw_path;
% 打印结果
disp(['最优迭代次数：' num2str(best_generation)])
disp(['最短长度:' num2str(1/fitness_best)])
disp 最佳路径:
tsp.print_path();
% 可视化结果
subplot(1,2,2);
plot(1./best_fitness)
save('./Best_Path_GA.mat','best_chromosome','-double');
